Spatial_data

datum = shape of the sphere CRS = is just a reference system area and distance measurements require a projection Projection is not always needed

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

library(leaflet)

Prepare the regions shapefile

ak_regions <-read_sf("data/shapefiles/data/shapefile_demo_data/ak_regions_simp.shp")
st_crs(ak_regions)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

class(ak_regions)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

ak_regions_3338 <- ak_regions %>% 
  st_transform(crs=3338)

plot(ak_regions_3338)

summary(ak_regions_3338)

##    region_id     region            mgmt_area          geometry 
##  Min.   : 1   Length:13          Min.   :1   MULTIPOLYGON :13  
##  1st Qu.: 4   Class :character   1st Qu.:2   epsg:3338    : 0  
##  Median : 7   Mode  :character   Median :3   +proj=aea ...: 0  
##  Mean   : 7                      Mean   :3                     
##  3rd Qu.:10                      3rd Qu.:4                     
##  Max.   :13                      Max.   :4

summary(ak_regions_3338)

##    region_id     region            mgmt_area          geometry 
##  Min.   : 1   Length:13          Min.   :1   MULTIPOLYGON :13  
##  1st Qu.: 4   Class :character   1st Qu.:2   epsg:3338    : 0  
##  Median : 7   Mode  :character   Median :3   +proj=aea ...: 0  
##  Mean   : 7                      Mean   :3                     
##  3rd Qu.:10                      3rd Qu.:4                     
##  Max.   :13                      Max.   :4

ak_regions_3338 %>% 
  select(mgmt_area)

## Simple feature collection with 13 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -2175328 ymin: 405653.9 xmax: 1579226 ymax: 2383770
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 13 x 2
##    mgmt_area                                                       geometry
##        <dbl>                                             <MULTIPOLYGON [m]>
##  1         3 (((-1156666 420855.1, -1159837 417990.3, -1161898 416944.4, -~
##  2         4 (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146, 5690~
##  3         3 (((-339688.6 973904.9, -339302 972297.3, -339229.2 971037.4, ~
##  4         3 (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 654303.1~
##  5         2 (((561012 1148301, 559393.7 1148169, 557797.7 1148492, 555907~
##  6         3 (((115112.5 983293, 113051.3 982825.9, 110801.3 983211.6, 108~
##  7         4 (((-678815.3 1819519, -677555.2 1820698, -675557.8 1821561, -~
##  8         4 (((-1030125 1281198, -1029858 1282333, -1028980 1284032, -102~
##  9         2 (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186, 36741~
## 10         4 (((-848357 1636692, -846510 1635203, -840513.7 1632225, -8358~
## 11         2 (((426007.1 1087250, 426562.5 1088591, 427711.6 1089991, 4299~
## 12         1 (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, 1296~
## 13         4 (((-375318 1473998, -373723.9 1473487, -373064.8 1473930, -37~

Prepare the population by region

pop <-read.csv("data/shapefiles/data/shapefile_demo_data/alaska_population.csv",stringsAsFactors = F)

class(pop)

## [1] "data.frame"

pop_4326 <- st_as_sf(pop,
                     coords =c("lng","lat"),
                     crs = 4326,
                     remove =F)# pass the crs that is currently is, not the one you want

pop_3338 <- st_transform(pop_4326,crs=3338)
#pop_3338<- st_coordinates() #add coordinate,
#add a colum with the ESPG

pop_joined <- st_join(pop_3338, ak_regions_3338, join=st_within)# we want to know which of the points is within each region
head(pop_joined)

## Simple feature collection with 6 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##   year     city      lat       lng population region_id           region
## 1 2015     Adak 51.88000 -176.6581        122         1 Aleutian Islands
## 2 2015   Akhiok 56.94556 -154.1703         84         6           Kodiak
## 3 2015 Akiachak 60.90944 -161.4314        562         8        Kuskokwim
## 4 2015    Akiak 60.91222 -161.2139        399         8        Kuskokwim
## 5 2015   Akutan 54.13556 -165.7731        899         1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153        777        13            Yukon
##   mgmt_area                   geometry
## 1         3  POINT (-1537925 472627.8)
## 2         3 POINT (-10340.71 770998.4)
## 3         4  POINT (-400885.5 1236460)
## 4         4  POINT (-389165.7 1235475)
## 5         3 POINT (-766425.7 526057.8)
## 6         4  POINT (-539724.9 1456223)

pop_region <- pop_joined %>% 
  as.data.frame() %>% 
  group_by(region) %>% 
  summarise(total_pop =sum(population))

head(pop_region)

## # A tibble: 6 x 2
##   region           total_pop
##   <chr>                <int>
## 1 Aleutian Islands      8840
## 2 Arctic                8419
## 3 Bristol Bay           6947
## 4 Chignik                311
## 5 Cook Inlet          408254
## 6 Copper River          2294

Now we need to get it back to the correct geometries.

pop_region_3338 <- left_join(ak_regions_3338, pop_region)

## Joining, by = "region"

# you could specify, which you would need if the colums did not have the same header
#pop_region_3338 <- left_join(ak_regions_3338, pop_region, by=c("region"="region")

#dissolve the boundaries 
pop_mgmt <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarise(total_pop=sum(total_pop))

plot(pop_mgmt["total_pop"])

# to keep internal boundaries 

pop_mgmt <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarise(total_pop=sum(total_pop), do_union=F)


plot(pop_mgmt["total_pop"])

pop_mgmt$ESPG = "3338"
?sf::tidyverse

## starting httpd help server ... done

Make maps in ggolot!

ggplot()+
  geom_sf(data = pop_region_3338, mapping = aes(fill=total_pop))+
  #geom_sf(data = rivers_3338, aes(size = StrOrder), color = "black") +
  geom_sf(data = pop_3338, aes(), size = .5) +
  scale_size(range = c(0.01, 0.2), guide = F) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma)

# Adding base maps in ggmap

3857 pseudomercador is the project on the tiels

pop_3857 <- pop_3338 %>% 
  st_transform(crs =3857)

# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}

bbox <- c(-170, 52, -130, 64)
ak_map <- get_stamenmap(bbox,zoom = 4)

## Source : http://tile.stamen.com/terrain/4/0/4.png

## Source : http://tile.stamen.com/terrain/4/1/4.png

## Source : http://tile.stamen.com/terrain/4/2/4.png

## Source : http://tile.stamen.com/terrain/4/0/5.png

## Source : http://tile.stamen.com/terrain/4/1/5.png

## Source : http://tile.stamen.com/terrain/4/2/5.png

ak_map_3857 <- ggmap_bbox_to_3857(ak_map)

ggmap(ak_map_3857)+
  geom_sf(data=pop_3857, aes(color=population),inherit.aes=F)+
  scale_color_continuous(low="khaki", high ="firebrick", label=comma)

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Leaflet

leaflet will project your files for you but you need to make sure you use somehting lik wgs 84 (4326) in needs to be a crs with lat long

# this would be a good function 
epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))

pop_region_4326 <- pop_region_3338 %>% 
  st_transform(crs=4326)

leaflet(options=leafletOptions(crs= epsg3338)) %>% 
  addPolygons(data=pop_region_4326,
              fillColor = "gray",
              weight =1)

#add legend scale 
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),# tells how to map it to the pallette created above
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft", #adds a legend to the map
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")
m

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1) %>% 
        addCircleMarkers(data = pop_4326,
                         lat = ~lat,
                         lng = ~lng,
                         radius = ~log(population/500), # arbitrary scaling
                         fillColor = "gray",
                         fillOpacity = 1,
                         weight = 0.25,
                         color = "black",
                         label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m